function [node, elem] = process_node_elem_2d(left, right, bottom, top, Nx, Ny, basis_type)

[X, Y] = meshgrid(linspace(left, right, Nx+1), linspace(bottom, top, Ny+1));
node = [X(:), Y(:)];

Np = size(node,1);
k = (1:Np-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
switch basis_type
    case {"P1", "P1b", "P2"}
        elem = [k, k+Ny+1, k+1;
            k+1, k+Ny+1, k+Ny+2];
        elem = sortrows(elem);

    case {"Q1", "Q1b", "Q2"}
        elem = [k, k+Ny+1, k+Ny+2, k+1];

    otherwise
        error("Invalid basis type");
end

switch basis_type
    case {"P1b", "Q1b"}
        new_node = zeros(size(elem,1),2);
        for k = 1:size(elem,1)
            V = [node(elem(k,:),:)];
            new_node(k,:) = polygon_centroid(V);
        end
        node = [node; new_node];
        new_elem_col = (Np+1:Np+size(elem,1))';
        elem = [elem, new_elem_col];

    case {"P2", "Q2"}
        [X, Y] = meshgrid(linspace(left, right, 2*Nx+1), linspace(bottom, top, 2*Ny+1));
        new_node = [X(:), Y(:)];
        cut = zeros(Ny+1,Nx+1);
        k = 1;
        for i = 0:2:2*Nx
            cut(:,k) = i*(2*Ny+1)+1:2:(i+1)*(2*Ny+1);
            k = k+1;
        end
        new_node(cut(:),:) = [];
        node = [node; new_node];
end

switch basis_type
    case "P2"
        t = Np+1;
        c = 1;
        for k = 1:2:size(elem,1)
            i = mod(c-1,Ny);
            new_elem_col(k:k+1,:) = [t+Ny+i, t+Ny+1+i, t;
                t+Ny+1+i, t+3*Ny+1, t+Ny+2+i];
            if mod(c,Ny) == 0
                t = t+2*Ny+2;
            else
                t = t+1;
            end
            c = c+1;
        end
        elem = [elem, new_elem_col];

    case "Q2"
        t = Np+1;
        c = 1;
        for k = 1:size(elem,1)
            i = mod(c-1,Ny);
            new_elem_col(k,:) = [t+Ny+i, t+3*Ny+1, t+Ny+2+i, t, t+Ny+1+i];
            if mod(c,Ny) == 0
                t = t+2*Ny+2;
            else
                t = t+1;
            end
            c = c+1;
        end
        elem = [elem, new_elem_col];
end

end